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The extragalactic sub-mm population: predictions for the 
SCUBA Half-Degree Extragalactic Survey (SHADES) 



o 
o 

(N 
X) 

00 

(N 
> 
(N 

in 
in 
oo 

o 
o 



o 



X 



Eelco van Kampen^'^, Will J. Percival^, Miller Crawford^, James S. Dunlop^, 
Susie E. Scott ^, Neil Bevis"^, Seb Oliver"^, Frazer Pearce^, Scott T. Kay"^, Enrique 
Gaztanaga^'^, David H. Hughes^, Itziar Aretxaga^ 

^ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ 
^ Institute for Astrophysics, University of Innsbruck, Technikerstr. 25, A-6020 Innsbruck, Austria 
^ Astronomy Centre, University of Sussex, Falmer, Brighton BNl 9QH 

* School of Physics & Astronomy, University of Nottingham, University Park, Nottingham, NG7 2RD 

^ Instituto Nacional de Astrofisica, Optica y Electronica, Apt. Postal 51 y 216, Puebla, Pue, Mexico 

^ Institut d'Estudis Espacials de Catalunya, Edifici Nexus, Gran Capita 2-4, desp. 201, 08034 Barcelona, Spain 

Accepted ... Received in original form ... 

ABSTRACT 

We present predictions for the angular correlation function and redshift distribution for 
SHADES, the SCUBA HAlf-Degree Extragalactic Survey, which will yield a sample 
of around 300 sub-mm sources in the 850 micron waveband in two separate fields. 
Complete and unbiased photometric redshift information on these sub-mm sources will 
be derived by combining the SCUBA data with i) deep radio imaging already obtained 
with the VLA, ii) guaranteed-time Spitzer data at mid-infrared wavelengths, and iii) 
far-infrared maps to be produced by BLAST, the Balloon-borne Large- Aperture Sub- 
millimeter Telescope. Predictions for the redshift distribution and clustering properties 
of the final anticipated SHADES sample have been computed for a wide variety of 
models, each constrained to fit the observed number counts. Since we are dealing with 
around 150 sources per field, we use the sky-averaged angular correlation function to 
produce a more robust fit of a power-law shape w{9) = {9/A)~^ to the model data. 
Comparing the predicted distributions of redshift and of the clustering amplitude A 
and slope 5, we find that models can be constrained from the combined SHADES data 
with the expected photometric redshift information. 
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1 INTRODUCTION 

The understanding of galaxy formation and evolution is 
making rapid progress because of a combination of a wealth 
of new observational data at a wide range of wavelengths and 
redshifts, and an increased understanding of which physical 
processes underlie the formation and evolution of galaxies. 
However, an important problem with most current galaxy 
formation models is that it is difficult to establish whether 
a set of model parameters that produces a good match to 
observations is unique. The main reason for this is that, for 
most models, there are degeneracies amongst the various 
free parameters. As a significant fraction of observational 
data used to constrain the model parameters are obtained 
from our local universe the 'uniqueness problem' can be re- 
solved by comparing model predictions and observations at 
high redshift, which in many respects is independent from a 
comparison at low redshift. 

One major difference between low and high redshifts 
is the frequency and intensity of major mergers: at high 



redshifts they occur far more often and are more intense as 
the participating galaxies are likely to be more gas-rich than 
their low-redshift counterparts. They are also expected to be 
dust-enshrouded at the peak of their burst of star formation, 
which means that they are most easily detected in the sub- 
mm or far-infrared wavebands. 

For this purpose, highly valuable observational data will 
be provided by the SCUBA HAIf-Degree Extragalactic Sur- 
vey (SHADES, see tittp://www. roe. ac.uk/ifa/ shades and 
Mortier et al. 2005 for details). This survey, which com- 
menced in December 2002, has been designed to cover 0.5 sq. 
degrees to a 3.5-a detection limit of Sgsoum ~ 8mJy, split 
between two 0.25 sq. degree fields. The two survey areas, 
the Lockman Hole East, and the Subaru-XMM Deep Field 
(SXDF), have been selected on the basis of low galactic con- 
fusion at sub-mm wavelengths, and the wealth of existing or 
anticipated supporting multi-frequency data from radio to 
X-ray wavelengths. 

In addition the SCUBA data will be combined with 
data from the VLA, the Spitzer telescope, and BLAST, 
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(see ittp: //chilel .physics .upenn. edu/blastpublic and 
Devlin 2001 for details), which will undertake a series of 
nested extragalactic surveys at 250, 350 and 500 pim. This 
experiment will significantly extend the wavelength range, 
sensitivity, and area of existing ground-based extragalactic 
sub- mm surveys (Hughes et al. 2002). 

It is anticipated that spectroscopic redshifts will ulti- 
mately be obtained for a substantial fraction of the SHADES 
sources (e.g. Chapman et al. 2003, 2004). However, the 
key point for the work presented here is that, even where 
optical/near- infrared spectroscopy is impossible, the long- 
wavelength data provided by the combined SCUBA -I- VLA 
-I- Spitzer -I- BLAST dataset will yield photometric redshifts 
for all sources with uncertainties o{ Sz < 0.5 (Aretxaga et 
al. 2004). This offers a unique powerful way of providing 
the complete and unbiased redshift and SED information 
required to measure the clustering properties of sub-mm 
sources, and the cosmic history of dust-enshrouded star for- 
mation that takes place in very massive star-bursts with 
inferred star-formation rates of order lOOOMoyr"^ (Scott 
et al. 2002). 

These massive star-bursts could be associated with the 
formation of the progenitors of massive ellipticals if sus- 
tained for a significant amount of time (up to 1 Gyr). How- 
ever, SCUBA sources could also be associated with bright, 
but short-lived bursts of intense star-formation occurring in 
more modest galaxies drawn from the high-redshift galaxy 
population already discovered at optical/UV wavelengths 
(Adelberger & Steidel 2000 and many others). If the bright 
SCUBA sources are indeed the progenitors of massive el- 
lipticals then they are likely to be more strongly clustered 
than when drawn from the population of less massive galax- 
ies. This is an inevitable result of gravitational collapse 
from Gaussian initial density fluctuations: the rare high- 
mass peaks are strongly biased with respect to the mass 
(Kaiser 1984). 

There is abundant evidence that this bias does occur 
at high redshift: the correlations of Lyman-break galaxies 
at a ~ 3 are almost identical to those of present-day field 
galaxies, even though the mass must be much more uniform 
at early times. Moreover, the correlations increase with UV 
luminosity (Giavalisco & Dickinson 2001), reaching scale- 
lengths of ro — 7.5 h~^Mpc, 1.5 times the present-day value. 
Daddi et al. (2000) find a trend of clustering with colour for 
EROs, reaching ro ~ 11 h^^Mpc for R — K > 5, which corre- 
sponds to fiuctuations in projected number density that are 
~unity on the scale of the SCUBAfield-of-view, falling to 
10% rms on 1-degree scales. For SHADES to detect the clus- 
tering properties of bright sub-mm sources over co-moving 
scales reaching >10Mpc, the survey needs to cover a sig- 
nificant fraction of a square degree. At the time of writing 
the survey is set to reach half a square degree within three 
years, and is making good progress towards achieving that 



The direct predecessor of SHADES was the 8-mJy sur- 
vey of Scott et al. (2002; see also Ivison et al. 2002). The cor- 
relation function for SCUBA sources derived from this sur- 
vey alone did not yield a significant detection of clustering, 
even though the large uncertainties meant it was still consis- 
tent with the strong clustering displayed by EROs. There are 
nevertheless good reasons for believing the SCUBA source 



population to be highly clustered, and some observational 
evidence for this is now found (Blain et al. 2004). In par- 
ticular, cross-correlation with X-ray sources (Almaini et al. 
2003) and Lyman-break galaxies (Webb et al. 2003) yield 
clearly-significant detections of clustering. 

Scott et al. (2005) have recently performed a combined 
clustering analysis on the three main existing blank-field 
SCUBA surveys (the 8-mJy survey, the CUDSS survey by 
Webb et al. 2003, and the Hawaii survey by Barger et al. 
1999) to determine whether the existing data are capable of 
revealing significant clustering within the sub-mm popula- 
tion alone. Even though this analysis is based on combining 
data from several small fields, it has yielded the first sig- 
nificant (5-0-) measurement of sub-mm source clustering on 
scales ~ 0.5—2 arcmin, of a strength that does indeed appear 
comparable to that found by Daddi et al. (2000) for EROs. 
Interestingly, if the integral constraint (see Section 4.2 for 
its definition) is varied as a free parameter, the inferred clus- 
tering in fact becomes stronger than that displayed by the 
ERO population. 

Finally, Blain et al. (2004) have found tentative evi- 
dence for a clustering length of (6.9 ± 2.1)/i^^ Mpc (comov- 
ing) for those submm galaxies for which they could obtain 
spectroscopic redshifts. As this was only possible for sources 
that are also detected at radio wavelengths, their sample is 
incomplete and likely to be biased. Furthermore, Adelberger 
(2004) argued that the method used is prone to systematic 
errors and is unnecessarily noisy. 

The paper is organized as follows: Section 2 provides 
an overview of SHADES and its main aims. Section 3 de- 
scribes the various models used to make predictions, which 
are compared to each other in Section 4. This section also 
presents the actual predictions for SHADES, and we discuss 
these in Section 5. 



2 SHADES: A WIDE-AREA SUB-MM SURVEY 
WITH REDSHIFT INFORMATION 

The science goals of SHADES are to help answer three fun- 
damental questions about galaxy formation: What is the 
cosmic history of massive dust- enshrouded star-formation 
activity? Are SCUBA sources the progenitors of present- 
day massive ellipticals? What fraction of SCUBA sources 
harbour a dust-obscured AGN? The aim of this paper is 
to review and compare the predictions of various existing 
models for the bright sub-mm population, and to consider 
how they can be best tested and constrained by the final 
SHADES dataset, and thus help answer the first two ques- 
tions. The third question it not addressed in this paper, as 
it involves a detailed analysis of the combined radio, mid- 
infrared and X-ray properties of the SHADES sources. 

An important property of SHADES is having mean- 
ingful redshift estimates, which provides vital information 
for estimating the bolometric luminosity of the sources, 
and hence the cosmic history of energy output from dust- 
enshrouded star-formation activity. Redshift information 
also holds the key to measuring the clustering properties 
of the sub-mm source population. Although precise spectro- 
scopic measurements of the redshift of a sample of SHADES 
sources will be possible if reliable radio/optical/IR counter- 
parts can be identified and readily foUowed-up with 10-m 
class optical telescopes (e.g. Chapman et al. 2003, 2004), 
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in practice one will not be able to derive this information 
for the majority of sources in the survey. However, combi- 
nation of the SHADES and BLAST data will allow the use 

of sub-mm photometric rcdshift techniques, yielding crude 
estimates (Sz < 0.5) for individual sources detected in both 
surveys. 

A Monte-Carlo based photometric-redshift technique 
has been designed by Hughes et al. (2002) and tested by 
Arctxaga ct al. (2003, 2004). Here sub-mm photometric in- 
formation is combined with prior information on the popu- 
lation, such as the number counts and the likely evolution 
of the luminosity function of dust-enshrouded galaxies, to 
weight the output redshifts provided by a large sample of 
template SEDs. These SEDs represent the wide range of 
temperatures, dust emissivities and luminosities found in 
nearby IR-bright galaxies. 

Even though the redshift distributions are relatively 
wide, the detailed information on the shape of the distri- 
bution, combined with a large number of sources, provides 
a powerful statistical measurement of population properties 
such as the parent redshift-distribution and the global star- 
formation rate. While, naively, these measurements might 
seem insufficiently crude, the combination of the redshift 
distributions of hundreds of sources can indeed measure the 
history of star formation of the galaxies detected in more 
than two sub-mm bands (Lfir > 2 x lO'^^ L©) with an ac- 
curacy of ~20% (Hughes ct al. 2002). 

While simulations show that photometric redshift esti- 
mates detected only from sub-mm data have errors of order 
0.5 (Hughes et al. 2002), it has been shown empirically that 
the inclusion of additional photometric information provided 
by detections or upper limits at 1.4 GHz (from the VLA), 
and at 170-70 /im (from the Spitzer Space Telescope), in- 
creases the accuracy of the photometric redshifts to ±0.3 
(Aretxaga et al. 2004). 

3 FOUR ALTERNATIVE MODELS OF THE 
EXTRAGALACTIC SUB-MM POPULATION 

Four different models for the clustering of SCUBA galaxies 
are presented: a 'simple merger' model, a 'hydrodynamical' 
model, a 'stable clustering' model and a 'phenomenological' 
model. Some of these are designed especially with SHADES 
in mind, while for other models the SCUBA predictions are 
just part of a range of predictions. The models also vary in 
the level of complexity, and in the underlying assumptions, 
including the choices for the cosmological parameters, even 
though differences in the latter axe minor compared to the 
fundamental differences between the models. 

The aim of this paper is not to perform a detailed com- 
parison between these models, or between models and data, 
but simply to present predictions for a diverse range of re- 
alistic models. The goal is to study the ability of SHADES 
to measure clustering, and establish its capability to distin- 
guish between models. 

3.1 A simple merger model 

The simple merger model is included in order to help de- 
termine the important processes at work in the creation of 
SCUBA galaxies. The underlying premise of this model is 
that 8-mJy SCUBA galaxies are formed by obscured star 



formation driven by the violent merger between two galaxy 
sized haloes. The emission is assumed to be above the 8-mJy 
detection threshold for a lifetime tufe after the galaxy haloes 

have merged. No direct link is made between the luminosity 
of the SCUBA galaxy and the properties of the merger ex- 
cept that a lower limit is placed on the final mass of haloes 
that contain a detectable 8-mJy SCUBA source. In other 
words, a Poisson sampling of massive halo mergers is as- 
sumed to form bright SCUBA galaxies. We have adopted a 
mass limit of lO'^^ M©, corresponding to 'radio-galaxy' mass 
haloes. 

Halo mergers were found in a 256'^ N-body simulation 
run within a co- moving (100/i~^Mpc)^ box using GADGET, 
a publicly available parallel tree code (Springel, Yoshida 
& White 2001). Cosmological parameters were assumed 
to have their concordance values (fim = 0.3, SIa = 0.7, 
h = 0.70 &i Us = 1), and the power spectrum normaliza- 
tion was set at as — 0.9. Outputs from the simulation were 
obtained at 434 epochs, separated approximately uniformly 
in time, and halos were found at each epoch using a stan- 
dard friends-of-friends routine with linking length b = 0.2. 
New halos were defined to be halos with > 50% of the con- 
stituent particles not having previously been recorded in a 
halo of equal, or greater mass. Of these, the halo was said 
to have been created by major merger if there were two pro- 
genitors at the previous time output that had mass between 
25% and 75% of the final mass. 

Obtaining the right number density of SCUBA sources 
is limited by the definition of merging used, the lifetime of 
emission above the detection threshold, and the proportion 
of mergers that result in SCUBA sources. We therefore sim- 
ply assume that all of the mergers, defined as above, result 
in a luminous SCUBA source, and allow tu^c to vary to give 
~300 sources in 0.5 deg^. Because the density of high-mass 
(> 10'^ Mo) mergers is low, obtaining the correct number 
density of SCUBA sources required a relatively long lifetime 
tiifg = 8 X 10* years. 

Mock SCUBA catalogues for a 0.5 deg^ survey were cal- 
culated by placing a (co-moving) light cone through an ar- 
ray of simulation boxes. This is done by selecting output 
time-steps such that corresponding redshifts are separated 
by a box length in co-moving coordinates. Boxes are re- 
fiected, rotated and translated randomly to reduce the ar- 
tificial correlation between neighbouring boxes inherent in 
using a single simulation, this necessarily reduces real cor- 
relations due to structures that would cross boxes. Mergers 
that occurred less than the model lifetime before the time 
corresponding to their luminosity distance were flagged as 
potential SCUBA sources, and their angular positions and 
redshifts were recorded in order to create mock catalogues. 

Obviously, while this model does predict both the 

spatial distribution and redshift space distribution of the 
SCUBA sources, it does not predict the luminosity func- 
tion. In fact, we note that following successful comparison 
between analytic theory and numerical simulations, both the 
redshift space distribution and the spatial distribution of 
SCUBA galaxies in this model could have been accurately 
estimated analytically (Percival, Miller & Peacock 2000; Per- 
cival et al. 2003). 
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3.2 A hydrodynamical model 

At the heart of the model is a simulation from Muanwong 
et al. (2002) that is an adaptive particle-particle, particle- 
mesh code incorporating smoothed particle hydrodynam- 
ics (SPH). The underlying code is HYDRA (Couchman, 
Thomeis & Pearce 1995) with the addition of a standard 
pair- wise artificial viscosity (Thacker et al. 2000). The cos- 
mological model is fim = 0.35, ^Ia = 0.65, h = 0.71, 
as = 0.9, Ob = 0.019/i~^. The simulation used here em- 
ploys a box of co-moving size (100/i~^Mpc)^ with 160^ dark 
matter particles and 160^ gas particles, and is evolved be- 
tween 50 > z > in approximately 2000 time-steps. The 
simulations have various components: non-interacting dark 
matter; gas; "star-like" (same as gas, but forming stars); 
"galaxy-fragments" , which are collisionless. Evolution of the 
various components is as follows: all particles evolve under 
gravity; gas can adiabatically heat and cool; gas can also 
radiatively cool; at p/p > 500, T < 12, OOOK gas particles 
become "star-like" (at this point all the mass is deemed to 
have been converted into stars). An aggregation of 13 or 
more close "star-like" particles become a "fragment". Frag- 
ments may accrete more star-like particles but do not merge. 

As a complete model of galaxy formation this simulation 
has a number of strengths and weaknesses. It provides a 
self consistent treatment of largo-scale-structure and galaxy 
evolution. However, the limited resolution and the arbitrary 
solution to "cooling catastrophe" necessitated by this, limit 
its validity. For the present purpose the full power of the 
simulation is not used; it serves as an ingredient to a more 
phenomonological model. 

The first step is to construct a "galaxy-fragment" light- 
cone in the usual way (as described in Section 3.1). With 
the SHADES sample area it is not necessary to use more 
than a single box transverse to the line-of-sight. 

The rcdshift distribution of the fragments in this cone 
{dN/dz)frag is measured using redshift bins of uniform 
width Az = 0.7, however there are many more fragments in 
each bin than would be detectable. Each fragment is treated 
as the possible location of a SCUBA source and is selected 
based upon its star-formation rate (SFR) , which is measured 
cis the mass per unit time of " star-like" particles accreted to 
this fragment (averaged over the last output time-step) . The 
required number of fragments with the greatest SFR are 
selected from each bin such that the redshift distribution 
matches a particular model: (dN / dz) scuba- For this pa- 
per, the analytical form of Baugh et al. (1996) was adopted 
with a median redshift of 2.3 and the normalization such 
to give 300 sources in the full 0.5 dcg^ sample size. Hence, 
the redshift distribution produced is not derived from the 
hydrodynamical simulations and this model merely makes 
a reasonable choice as to which fragments SHADES will in- 
clude, and for these, encodes the positional information from 
the simulations. 

3.3 A stable clustering model 

This is the model of Hughes & Gaztanaga (2000), in which 
a single output from a N-body simulation that fits well 
the local spatial correlation function as measured by APM 
(Gaztaiiaga & Baugh 1998) is used to generate a population 
of SCUBA galaxies in a lightcone. This corresponds to as- 
suming stable clustering, i.e. a constant spatial correlation 



function in co-moving space. Fixing the spatial correlation 
function does not imply that we also fix the angular correlar 
tion, as that depends on lightcone geometry, luminosity evo- 
lution, and the rcdshift selection function. The prescription 
for galaxy formation corresponds to the assumption that the 
probability for finding a galaxy somewhere within the light- 
cone is simply proportional to the local dark matter density, 
with the total number of galaxies normalized to the surface- 
density required for a given flux limit. Although the redshift 
distribution is (exponentially) cut-off beyond 2 = 6, this 
model contains the highest redshift SCUBA galaxies of all 
models considered in this paper. 

This model has been used in the photometric redshift 
estimation technique of Arctxaga et al. (2003), to constrain 
sample size and depth given the correlation length, and to 
test correlation function measurements from surveys with 
relatively small sky coverage (Gaztanaga & Hughes 2001). 

3.4 A phenomenological model 

The phenomenological galaxy formation model of van Kam- 
pen, Rimes & Peacock (2005), a revised version of the 
model of van Kampen, Jimenez & Peacock (1999), is semi- 
numerical, in the sense that the merging history of galaxy 
haloes is taken directly from N-body simulations that in- 
clude special techniques to prevent galaxy-scale haloes un- 
dergoing 'overmcrging' owing to inadequate numerical reso- 
lution. When haloes merge, a criterion based on dynamical 
friction is used to decide how many galaxies exist in the 
newly merged halo. The most massive of those galaxies be- 
comes the single central galaxy to which gas can cool, while 
the others become its satellites. 

When a halo first forms, it is assumed to have an 
isothermal-sphere density profile. A fraction fib /Om of this 
is in the form of gas at the virial temperature, which can 
cool to form stars within a single galaxy at the centre of the 
halo. Application of the standard radiative cooling curve 
shows the rate at which this hot gas cools and is able to 
form stars. Energy output from supernovae reheats some of 
the cooled gas back to the hot phase. When haloes merge, 
all hot gas is stripped and ends up in the new halo. Thus, 
each halo maintains an internal account of the amounts of 
gas being transferred between the two phases, and consumed 
by the formation of stars. 

The model includes two modes of star formation: qui- 
escent star formation in disks, and star-bursts during major 
merger events. Having formed stars, in order to predict the 
appearance of the resulting galaxy it is necessary to assume 
an IMF, which is generally taken to be Salpeter's, and to 
have a spectral synthesis code, for which we use the spec- 
tral models of Bruzual & Chariot (1993). The evolution of 
the metals is followed, because the cooling of the hot gas 
depends on metal content, and a stellar population of high 
metallicity will be much redder than a low metallicity one 
of the same age. It is taken as established that the popula- 
tion of brown dwarfs makes a negligible contribution to the 
total stellar mass density, and the model does not allow an 
adjustable M/L ratio for the stellar population. The cos- 
mological model adopted is Qm = 0.3, Oa = 0.7, h = 0.7, 
CTg — 0.93, Ob ~ 0.02/i~^. The 850/im fiux is assumed to be 
proportional to the star formation rate (with 8mJy corre- 
sponding to 1000 Me/yr, as found by Scott et al. (2002), 
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Figure 1. A simulated distribution of 8-mJy sources for tlic plic- 
nomenological model described in Section 3.4, with no rcdsiiift 
selection. The size of the field is equal to each of the two fields 
that will make up SHADES: a quarter square degree. The diam- 
eter of each dot is proportional to the sub-mm magnitude of the 
source it represents. 

with a random term of order 50 per cent added or subtracted 
to mimic the uncertainty in dust temperature, grain sizes, 
and other properties that are not yet included in the mod- 
elling. 

The model used in this paper has a mixture of bursting 
and quiescent star formation, with most of the recent star 
formation occurring in discs, following the Schmidt law with 
a threshold according to the Kennicutt criterion, and most 
of the high-redshift star formation resulting from merger- 
driven star-bursts. The model is similar in philosophy to that 
of Hatton et al. (2003) and Baugh et al. (2004), but with 
many differences in the details and parameters adopted. 

3.5 Model comparison 

In this section we compare the models in order to get a 
qualitative description of where the differences lie. 

What drives the flux in each of the models? In the sim- 
ple merger model, an actual flux is not calculated; rather, a 
merger mass threshold is related to a flux threshold. The hy- 
drodynamical and stable clustering models have some of the 
properties of the galaxies fixed, but not the flux, which is as- 
signed statistically. The phenomenological model is the only 
one that generates a flux from the actual physical properties 
of the galax;ies, taking into account the approximations and 
assumptions made. 

What drives the clustering signal in each of the mod- 
els? In the stable clustering model the spatial correlation 
function is fixed in comoving space, which means that the 
angular correlation function is built up along the lightcone 
in a way that depends on the selection of galaxies as a func- 
tion of redshift. The simple merger, hydrodynamical and 
phenomenological models produce galaxies first, and built 
up the angular correlation function along the fine of sight. 
In the simple merger a one-to-one correspondence between 



galaxies and haloes is assumed, whereas the other two mod- 
els have more complex relations between mass and light. 

What determines the redshift distribution in each of 
the models? In the hydrodynamical model it is simply taken 
from Baugh et al. (1996), whereas for the other models it 
is actually an outcome of the models, although in the sta- 
ble clustering model an exponential cut-off at z « 6 is ap- 
plied. The main determining factor for the simple merger 
and phenomenological models for the redshift distribution 
is the merger rate as a function of redshift. For the simple 
merger model this is obvious, but for the phenomenologi- 
cal model this follows from the dominant contribution of 
merger-driven starbursts to the sub-mm flux. 

Besides the models used in this paper, other models ex- 
ist in the literature, which are similar in philosophy to those 
included here, but still produce different predictions for the 
sub-mm population. Models similar to the phenomenologi- 
cal model are those of Hatton et al. (2003) and Baugh et 
al. (2004), which differ mainly in the details of the physics 
implemented, and the choice of parameters. 



4 MODEL PREDICTIONS COMPARED 

In Fig. 1 we show a complete simulation of one of the 
two 850 fim datasets that will comprise SHADES, produced 
using the phenomenological model of van Kampen et al. 
(2004). A simple square geometry was chosen, although the 
actual survey geometry of each of the two SHADES fields 
could be of a somewhat different shape. All sources with 
fluxes larger than 8 mjy are shown, where the symbol size 
is proportional to the logarithm of the flux, i.e. a sub-mm 
magnitude. 

We now consider simulated SHADES datasets as pre- 
dicted from the various alternative models of the sub-mm 

source population presented in Section 3. We do not con- 
sider the effects of noise and sidelobes, but we do take into 
account the effects of the 15 arcsec SCUBA beam by merg- 
ing into single sub-mm sources anything closer together than 
7 arcsec. This reflects the resolution expected for the final 
source extraction from the SHADES images. After the re- 
moval of close pairs we assume the model sub-mm sources 
to refiect the final SHADES source list. 

The final survey is expected to contain around 300 
sources, i.e. 150 sources per field. We produce, for each 
model, 50 realizations of individual fields, i.e. 25 mock 
SHADES datasets of 300 sources each. 

4.1 Predictions for the redshift distribution 

For the four models, wc show in Fig. 2 the redshift distri- 
butions expected after smoothing with a Gaussian filter of 
radius 0.4, which reflects, very crudely, the resolution achiev- 
able with photometric redshifts. The distributions are ob- 
tained by averaging over all realisations for each model, and 
are normalized to the total source count of 300. 

Even after relatively heavy smoothing, we see that the 
redshift distributions are rather different, and are clearly dis- 
tinguishable from each other. This means that even crude 
but complete redshift information will be of enormous bene- 
fit in differentiating between and constraining models. Obvi- 
ously, obtaining more accurate redshifts should help to tune 
the models that survive this first test even further. 
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Figure 2. Redshift distributions for all models with 300 sources, 

averaged over 25 mock SHADES datasets for each model. The 
model distributions were convolved with a Gaussian of radius 0.4, 
roughly mimicking the statistical uncertainty in the photometric 
redshift errors. 



For the purpose of comparing clustering properties be- 
tween the models, note from Fig. 2 that the redshift range 
2 < a; < 3 is the only range where all models have a rea^ 
sonable number of sources to attempt a correlation func- 
tion analysis. The differences in the distribution stem from 
the different assumptions for each of the models: the sim- 
ple merger model assumes that only high-mass mergers can 
form SCUBA sources, which, in a hierarchical structure for- 
mation scenarion, necessarily places them at lower redshifts 
as compared to the other models, while the simple, unbiased 
galaxy formation prescription of the stable clustering model 
places SCUBA sources at relatively high redshifts. 

Current spectroscopic redshift measurements for sub- 
mm selected galaxies are incomplete and only available for 
small samples, so any redshift distribution derived from such 
measurements is tentative. Chapman et al. (2004) claim that 
a Gaussian distribution with z = 2.4 and ffz = 0.65 fits their 
available data well, but the incompleteness of their dataset 
imposed through their radio selection hinder a comparison 
to models. 



4.2 Clustering measures 

The estimated redshifts have a predicted accuracy of (52 ~ 
±0.4, which means that we cannot directly measure the 3D 
spatial correlation function ^(r). However, we do not have to 
restrict ourselves to measuring angular clustering, as photo- 
metric redshifts can be used to boost the angular clustering 
signal-to-noise by splitting the sample in redshift bins, or by 
only considering pairs of galaxies that lie at similar redshifts. 
Even so, the measured correlation function will be noisy, so 
we use integrals of this function, as considered in the early 
days of optical gala^xy surveys when total source counts were 
much lower than today (e.g. Davis & Peebles 1983). 



4-2.1 Estimating the angular correlation function 

The method for modelling the clustering of sources proceeds 
as follows. From the data, the Landy & Szalay (1993) esti- 
mator 1 + wls = 1 + {DD - 2DR + RR)/RR is calculated, 
where DD, DR and RR are the (normalized) galaxy-galaxy, 
galaxy-random and random-random pair counts at separa- 
tion 6, calculated from the galaxy sample and a large random 
catalogue containing 10000 points that Poisson samples the 
survey region. This estimator is then fitted by its expected 
value 



l + {wLs) = [l-hM)(e)]/(l-Fwn) 



(1) 



where wn is the integral of the model two-point correlation 
function over the sampling geometry: 



Jn 



(2) 



The function Gp{6) is the probability density function of 
finding two randomly placed points in the survey at a dis- 
tance 6. This "integral constraint" corrects for the effect of 
not knowing the true density of objects (Groth & Peebles 
1977; Landy & Szalay 1993) and stops the recovered correla- 
tion function being biased to low values compared with the 
true function. Note that cq. (1) implies that the true cor- 
relation function is biased low by a factor 1 + , whereas 
often this is approximated by w{d) = WLS—wn, i.e. ignoring 
the term WLSWn- 

Wc also introduce an alternative to the standard angu- 
lar correlation function that takes redshift information into 
account in an unorthodox way. In the counting of DD pairs, 
we just consider those pairs that have a redshift separa- 
tion of at most 0.4, whereas the DR and RR counts are 
still obtained for all galaxy pairs. This is equivalent to re- 
moving distant pairs that are expected to be unclustered 
from an analysis of the angular correlation function of all 
of the objects in the survey. It is clear that this approach 
must increase the signal-to-noise of the recovered correlation 
function. 



4-2.2 The sky-averaged angular correlation function 

For the relatively small number of sources being detected in 
SHADES, we measure an integral of the Landy & Szalay es- 
timator. Such an approach has previously been used to anal- 
yse clustering within early galaxy redshift surveys (eg. the 
CfA survey; Davis & Peebles 1983). The statistic that was 
often obtained in these analyses was the integrated quantity 
Ja, defined as 

Mr) = f i{y)y'dy . (3) 
J a 

The dimensionless analogue of J3 is called the volume- 
averaged correlation function: 



e(0 



Mr] 



(4) 



This measures the fluctuation power up to the scale r, and 
is therefore a useful measure for a survey that is limited 

in object numbers. For reference, ^(10ft~^Mpc) = 0.83 was 
found for the optical CfA survey (Davis & Peebles 1983; no 
error given). 
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As we cannot measure the spatial clustering function 
^(r), as the redshift determinations axe very uncertain, we 
use the two-dimensional version of f , the sky-averaged an- 
gular correlation function 

2 r" 

w{e) = J wi<P)<pdct> , (5) 

where w{6) is the angular correlation function, which is the 
projection of ^(r) along the line-of-sight. Our estimator of 
this statistic was calculated by numerically integrating the 
angular correlation function (calculated using the Landy & 
Szalay 1993 estimator), in the form of logarithmically binned 
estimates Wi, up to the angle 9i using eq. (5): 

Wi = -^^WjOjA , (6) 

where A is the logarithmic binsize. The errors on Wi are 

obtained by propagating the errors on Wi through this sum- 
mation. This estimate for the true sky-averaged angular cor- 
relation function is also biased, and has its own integral 
constraint, similar to the one for w{d). For a power-law cor- 
relation function ■Wpi{6) — {d/A)~^ 

Thus, the sky-averaged correlation function is also a power- 
law, with the same slope and different amplitude (except 
for 5 = 2), and the integral constraint wq is scaled by the 
same factor 2/(2 — 5) with respect to wn (see eq. 2). We can 
therefore fit power- law models to either 'w{0) or w(0), and 
constrain the same parameters. 

4-2.3 Fitting to the model correlation functions 

Traditionally, minimization is used to fit a power-law 
function to the estimated correlation function. This type of 
minimization is, strictly speaking, only valid for binned data 
that is uncorrelated and has Gaussian errors. Our estimates 
of w{ff) and ■w{ff) are correlated for different maximum separ 
rations, and have errors that are non-Gaussian (it is easy to 
sec that, in the absence of clustering, the errors arc strictly 
Poisson, as shown by Landy & Szalay (1993) for the estima- 
tor of w(0)). In order to constrain models of the correlation 
function using our binned estimates, we should therefore 
perform a full likelihood calculation taking into account the 
potentially complex shape of the likelihood. Both w{6) and 
'w{6) obviously contain the same information, and should 
therefore result in the same likelihood surface for a given 
model. 

The reason for fitting w{9) rather than wiO) lies in 
the approximations that are made in estimating the like- 
lihood. For the sky-averaged angular correlation function, 
the bins are dependent on more pairs of galaxies than the 
corresponding direct estimator of the correlation function, 
so that the sky-averaged statistic will have a distribution 
closer to a Gaussian form. Switching from a direct estimate 
to a sky-averaged estimate increases correlations between 
data points. However, this can be taken into account prop- 
erly in the fitting procedure by calculating the full covari- 
ance matrix for the binned data, which is then diagonal- 
ized by a unitary transformation to produce an alternative 



X^-statistic (e.g. Fisher et al. 1994). This statistic is subse- 
quently employed to fit models to the data. 

A single parameter fit to the correlation function is of- 
ten adopted assuming that ■w{e) = {9/A)-°-^ (e.g. Roche et 
al. 1993; Daddi et al. 2000). It is not at all clear what the 
slope for high-redshift sub-inin sources is going to be, but it 
is easier to fit a one-parameter function for a small number of 
galaxies. In the following we consider both a one-parameter 
model with constrained power-law slope S = 0.8, and a two- 
parameter fit for the generic power-law w{9) = (9/A)~^ 
to both w{9) and w{9). We employ non-linear x^-fitting 
for both functions, using the Levenberg-Maxquaxdt method 
(Press et al. 1988), which allows us to easily take into ac- 
count the multiplicative integral constraint. For each fit the 
probability Q is calculated using the incomplete gamma 
function, and any fits with Q < 0.1 are discarded. 

4-2.4 Examples of mock angular correlation functions 

In Figs. 3, 4 and 5 we show examples of w{9) and w{9) for 
realisations of all models, for the case where we have no 
redshift information (Fig. 3) and for the case where we have 
(Fig. 4 and 5). The same realisations have been used for 
all three figures. Round symbols show the angular correla- 
tion function w{9), where open symbols indicate negative 
values, and error bars indicate Poisson errors. The best fit 
power-law is shown as a solid lino, if a fit was possible. If 
not, the line is simply omitted. The sky-averaged angular 
correlation function w{9) is shown using stars, and the best- 
fitting power-law is shown as a dashed line (again, if a fit 
was successful). 

Before continuing, please note that those examples arc 
by no means meant to be representative, they are merely 
shown to demonstrate the difference between using w{9) and 
w{6) for fitting, and to show the effect of the additional 
redshift information. The examples should not be used to 
compare the differences in clustering strength between the 
models; this will be covered in the next section. 

First focusing on Fig. 3, we see that the correlation 
functions for the complete linc-of-sight are noisy, and for 
two of the models a fit to w(0) fails completely. The figure 
demonstrates the use of sky-averaging, as w(9), plotted using 
stars, is better behaved. This is perhaps best illustrated in 
the first panel, but also in the third panel, where it has the 
direct consequence that a fit to w{9) is possible where a fit 
to w{9) failed (third panel of Fig. 3). In general, for most 
realisations a simple fit to ■w{9) turns out to be easier 
than a fit to the angular correlation function itself. This 
is very helpful for our purpose of comparing models, as it 
is important that fitting is possible for a large number of 
realisations, which needs to be largely automatic. 

A stronger clustering signal is expected for sources se- 
lected in a redshift range, as the signal is less polluted by 
uncorrelated sources at very different redshifts. Indeed, Fig. 
4 shows that for all realisations the selection of sources in 
the redshift interval 2 < z < 3 boosts the clustering signal. 
Even though the errors are larger due to the small number of 
sources, the stronger signal means that a fit is possible in all 
cases shown, even for w(9) itself, although the sky-averaged 
correlation function is to be preferred nevertheless. However, 
a fraction of realizations still produce unacceptable fits. 

An alternative to redshift intervals is to only count 
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Figure 3. Bias-corrected angular correlation function w{0) 
(round symbols, open meaning negative) and its sky-averaged 
counterpart w(0) (stars, negative values not plotted) for a sin- 
gle realisation of each model, with best fit power-law functions 
overplotted for both, and no redshift selection. See text for full 
details. 



Figure 4. Same as Fig. 3, but with sources selected to be in the 
redshift range 2 < z < 3, which is the range where all models 

have a reasonable amount of sources in this redshift bin (sec Fig. 
2). All realisations are exactly the same as in Fig. 3, in order to 
demonstrate what redshift availability can achieve. 
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Figure 5. Same as Fig. 3, but with sources selected to be in 
redshift pairs with dz < 0.4. All realisations are exactly the same 
as in Figs. 3 and 4. 

galaxies that are paired in redshift space, e.g. have \zi~Zj \ < 
0.4, as described in Section 4.2.1. The result for the same 
realizations as used for Figs. 3 and 4 is shown in Fig. 5, 
and again much better results are obtained as compared to 
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Figure 6a. Scatter plot for the two fitting parameters of the 
sky-averaged angular correlation function, A^i^y and S^j^y, for 50 
fields of 150 sources. Only fits of sufficient quality, i.e. those with a 
probability larger than 0.1, are included (see text for details). 



the estimates without any redshift information (Fig. 3). The 
binned data looks cleaner than those for the redshift inter- 
vals (Fig. 4), which is due to more galaxies being used in the 
DD counts. The fitting is therefore somewhat more reliable, 
as demonstrated by the small difference between fits to w(9) 
and id{9). 

4-2.5 Distribution over fitting parameters 

So far we have considered single realisations for each model, 
which should really be treated as examples of how the final 
SHADES dataset might appear. In order to be able to make 
a quantitative comparison possible between the actual final 
SHADES dataset and all models, we need to find the probar 
bility distribution over the fitting parameters for each model 
given the survey constraints (area, flux limit, etc.). 

We therefore produced 50 realizations for each model, 
and fitted a power-law correlation function to all of these. 
The resulting amplitudes ^sky and slopes Ssky of the sky- 
averaged correlation function w{8) are shown as scatter plots 
in Fig. 6a, which shows the case where no redshift informa- 
tion is available, and in Figs. 6b and 6c, where we are able 
to split up the sample in redshift intervals (three of these 
are shown), or select pairs of galaxies in redshift space. 

In the case of no redshift information, a significant frac- 
tion of the mock fields do not produce a correlation function 
that can be fitted by a power-law, and the number of esti- 
mates in Fig. 6a is therefore less than 50 for each model. 
However, for each model still more than half the realiza- 
tion allow a good fit so, crudely speaking, one would expect 
that at least one of the two observed SHADES fields should 
produce a good fit. All models spread out over a relatively 
large region of parameter space, and seem to overlap with 
each other for most of that region. This merely reflects the 
fact that whatever intrinsic correlation exists in the under- 
lying sub-mm population is weakened by projection, which 
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Figure 6b. As Fig. 6a, but for tliree different redshift bins. Again 
only parameters from 'good' fits are shown, i.e. those with a 
probability larger than 0.1. 



Figure 6c. As Fig. 6a, but for close redshift pairs with 5z < 0.4. 
Again only parameters from 'good' fits are shown, i.e. those fits 
that have a probability larger than 0.1. 



produces this large spread in fitting parameters and the sig- 
nificant overlap between the models. Only the high-mass 
merger model shows a larger clustering amplitude overall, 
which also results in a larger fraction of the realizations pro- 
ducing good fits, and a smaller spread in the slope 6. 

Let us now use the redshift information to split up the 
mock samples into redshift intervals, which should show 
stronger clustering. Three intervals arc shown in Fig. 6b, 
where the 2 < z < 3 case is the most relevant one as it 
has the most sources for all models. The other two intervals 
both have at least one model where the number of sources is 
significantly lacking, which means that only the remaining 
models can reasonably be compared. The first thing to no- 
tice in all three panels of Fig. 6b is that the clouds of fitted 
parameter pairs start to separate out somewhat, reducing 
the overlap between the models. 

Various interesting effects can bo soon for the different 
models. The stable clustering model (open squares) shows 
fairly strong clustering in the 1 < z < 2 redshift interval, but 
because of the low number of sources in this redshift range 
(see Fig. 2) , few of the realisations actually produce a good 
fit. The high mass merger model (open diamonds) shows rel- 
atively strong clustering for 1 < « < 3, but of course lack 
numbers in the highest redshift interval. The phenomenolog- 
ical model (crosses) shows strongest clustering in that same 
3 < z < 4 interval, but this is also somewhat troubled by 
low source counts. It also shows the weakest clustering for 
the lowest redshift interval. The clustering strength of the 
hydrodynamical model (open triangles) is virtually indepen- 
dent of redshift. 

In Fig. 6c we show the results for the close pairs, i.e. 
galaxies with \Sz\ < 0.4. The points scatter in a similar 
fashion to the 2 < z < 3 interval, with some differences, and 
the separating of model point clouds is comparable between 
models, except for the simple merger model, which can quite 
clearly be distinguished. This diagram should thus be a good 
test for high-mass merging versus the other models consid- 
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Figure 7. Distribution of clustering amplitude Agj^y over fifty re- 
alisation for each model, for the redshift bin 2 < z < 3, where each 
model has a sufficient number of sources available (top panel) , and 
for close redshift pairs with Sz < 0.4 (bottom panel). 



Figure 8. Distribution of the angular correlation measure «)(1') 
over fifty realisation for each model, again for the redshift bin 
2 < z < 3 (top panel), and the close redshift pairs with Sz < 0.4 
(bottom panel). 



ered in this paper. 

If we concentrate on the 2 < z < 3 interval, which al- 
lows the cleanest comparison between the four models con- 
sidered here, and on the close pairs, we do see that the clouds 
of points are overlapping significantly, but the distributions 
are elongated somewhat along the 6 axis, and have different 
mean clustering amplitudes A. Also, the mean of the distri- 
bution is near 5 = 0.8, observationally found for a range of 
galaxy types. 

This leads us to finally consider the traditional one- 
parameter fit to the data, assuming S — 0.8. This produces 
a single clustering amplitude A^ky for each mock field, and 
a distribution over ^sky for each model. These distributions 
are plotted in Fig. 7, for the redshift interval 2 < z < 3 (top 
panel), and for the close pairs. Interestingly, for this redshift 
interval, all distributions are different, and although there is 
significant overlap, the final SHADES dataset will distin- 
guish between the models, especially in combination with 
the different redshift distributions (sec Fig. 2). For the close 
pairs, the result of Fig. 6c is made more apparent, in that 
the high-mass merger model is clearly different from the rest 



of the models, which show almost identical distributions. 

Another measure of clustering is the sky-averaged angu- 
lar correlation function as it was originally intended: just as 
a measurement. Therefore we also plot, in Fig. 8, the distri- 
bution over a particular w{Oi), which we choose to be w(l'), 
i.e. the sky-averaged correlation function within 1 axcmin, 
for the same redshift selected data as used for Fig. 7. The 
result is a similar, although the distributions overlap more 
than those for ^^ky (as seen in Fig. 7). However, the major 
advantage with the measure w(l'), or one at a different an- 
gle, is that we do not need to assume a model for the form 
of the correlation function. 

5 CONCLUSIONS 

One of the primary science drivers for the SHADES project 
is to place strong constraints on galaxy formation models by 
observations of luminous sub-mm galaxies in the high red- 
shift Universe. In order to achieve this, it is worth consider- 
ing how models are best constrained by the data, and exam- 
ine the range of possible predictions. With this aim in mind, 
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we have presented 4 different models of the sub-mm galaxy 
population, selected to be widely varying in concept, without 
worrying about every aspect of each model. We avoid any 
strong assumption about the nature or redshift distribution 
of the sub-mm population. Given the uncertainties in what 
is known about the sub-mm population, we want to keep 
open a range of models, even those that can be questioned 
in some of their aspects. 

For each model, the redshift distribution and cluster- 
ing properties of the sub-mm population were predicted 
for SHADES, and 50 realisations were produced, each 
comprising around 150 sources. Thus, 25 mock SHADES 
datasets have been produced for each model. These simu- 
lated SHADES catalogues were used to investigate the abil- 
ity of the clustering statistics of the final dataset to constrain 
the various models collected here. Direct and sky-averaged 
estimators of the correlation function have been considered 
and their relative merits discussed. We have argued that 
power-law fits are best performed on the sky-averaged an- 
gular correlation function, and that a relatively good fit is 
possible in most cases for this measure. The full covariance 
matrix was calculated and diagonalized resulting in a 
statistic that was used to fit the power-law model to the data 
using the Levenberg-Marquardt method as implemented by 
Press et al. (1988). 

All models predict sufficiently strong clustering, so that 
we expect to detect clustering within the SCUBA popula- 
tion when SHADES is complete. Although cosmic variance 
remains a concern, we can certainly quantify the probability 
of a given model to produce the observed dataset, and this 
is expected to reject some of the models included in this pa- 
per. However, the aim of this paper is not to constrain the 
models; this will be done when the full SHADES dataset 
is available. In fact, we simply have observed basic trends 
between models by reducing each measured correlation func- 
tion to two power-law parameters: its slope and amplitude 
respectively. 

The observed redshift distribution will provide an com- 
plementary strong test of models, even with relatively coarse 
photometric redshift information. In fact, the combination of 
clustering and redshift data offers the best discriminator be- 
tween the different models that we have considered: models 
with similar redshift distributions have different clustering 
strengths, while models with similar clustering properties 
have different redshift distributions. 

Recently, Blain et al. (2004) considered the clustering 
of sub-mm sources in a number of relatively small fields for 
data with follow-up spectroscopic redshifts. An approach 
was adopted that selected galaxy pairs based only on radial 
positions, and did not use any angular information. Pairs 
of sub-mm sources with Av = 1200 km s~^, equivalent to 
separations of order 5 Mpc (comoving) a,t z — 2.5, were 
counted, and compared with the integrated model ^(r). For 
a larger survey where there is significant angular informa- 
tion, this method is not optimal, so for SHADES we prefer 
to fully exploit the angular information, with the restricted 
photometric redshift information as a subsample selection 
tool. 

The Blain et al. (2004) method is equivalent to the 
standard method of calculating ^(r) by pair counting (using 
DD/RR-1), but performing this in single bin with r < -Rmax. 
While angular clustering measurements ignore radial infor- 



mation, this "radial" method instead ignores angular in- 
formation in a similar way. The method should therefore 
include the equivalent of an integral constraint. Also, the 
quoted errors are derived from Poisson errors (on the pair 
counts), whereas the errors are expected to be larger than 
Poissonian for this estimator (Landy & Szalay 1993). Fur- 
ther limitations of the method employed by Blain et al. 
(2004) are discussed by Adelberger (2004). Given these is- 
sues, it is probably too early to use the Blain et al. clustering 
results to distinguish between models. 

The primary conclusion from our analysis is that, with 
the area coverage (0.5 square degrees) and the expected 
number of sources (200-400), and particularly with the 
expected photometric redshift information (Az ~ ±0.3), 
SHADES is capable of distinguishing between widely vary- 
ing scenarios for the production of the bright sub-mm pop- 
ulation. 
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